地震数据采集观测系统设计的重要目标之一是平衡成本和控制数据假频。即使是近年来发展起来的“两宽一高”地震数据采集, 对于中高频数据空间采样仍旧不够, 假频现象仍然不可忽视^[1]。假频的存在给地震数据去噪、偏移和反演都带来了很大困难。在偏移成像中, 地球物理学家分别对Kirchhoff积分偏移、单程波偏移、逆时偏移和高斯束偏移等成像技术的反假频进行了系列研究, 改善了相关算法的成像效果^[2-7]。假频通常划分为数据空间假频、算子假频和成像条件假频^[4]。不同的成像方法面对的假频也不相同, 如Kirchhoff积分偏移通常要处理算子假频^[5-6]、炮域单程波方程偏移和逆时偏移通常要处理数据假频和成像条件假频等^{[2, 4]}。当前, 几乎所有的工业偏移成像方法都会对假频进行处理^[3]。

双平方根方程偏移由于其计算效率和精度较高、偏移孔径无需人为选择、偏移角度道集容易输出等优点, 成为一种重要的成像方法^[8-10]。近年来, 业界研究人员持续研究该成像方法的各向异性介质和吸收衰减介质的推广、角度道集生成、传播精度改进以及实际数据的应用^[11-16]。对于双平方根方程偏移反假频研究, ZHANG等^[3]指出需要以半道间距成像来避免假频。在偏移实践中, 通常我们在中点-半偏移距域进行外推及成像^[8-12], 中点方向采样间隔为道距的一半, 看上去这种双平方根偏移不受假频影响。但是, 在实际偏移应用中, 需要同时满足中点-半偏移距两个维度(三维地震数据为四个维度, 分别为Inline方向和Crossline方向中点、Inline和Crossline方向半偏移距)均没有假频。ZHANG等^[3]仅指出以半道距成像还远远不够。对于双平方根方程偏移, 满足中点方向以半道间距采样情况下没有成像条件假频, 但对半偏移距方向的数据假频并没有展开论述, 这才是产生双平方根方程偏移假频的关键因素, 只保证波场以半道间距间隔传播和成像无法完全消除该算法的假频。对于当前三维实际数据, 炮线间距和检波线间距较大, 面元间距通常选择炮点和检波点间距的一半, 偏移距间距为炮线和检波线间距的两倍, 就造成一个面元中Inline方向和Crossline方向偏移距间距很大, 导致在中点-半偏移距域数据中假频现象非常严重。如果半偏移距方向的假频处理不到位, 双平方根方程波场延拓会引入大量噪声, 严重影响成像质量。这也是该方法多年来没有被工业界大规模应用的重要原因这一。因此, 开展双平方根方程偏移反假频的研究及消除十分必要。当前, 针对双平方根方程偏移的假频研究除了在文献[3]中简单提及, 并没有引起地球物理业界的广泛关注, 针对性的假频研究文献还较为少见。

本文首先简单回顾了双平方根方程偏移成像算法的基本原理, 然后从采样定理角度分析了成像过程中假频的来源, 给出了解决中点-半偏移距域数据空间假频问题的三角平滑滤波方法, 指出了该方法相较数据插值加密的优点。最后, 利用Marmousi模型数据验证了本文方法的有效性, 并对比了与常规方法计算需求的差异, 进一步利用实际数据对本文方法进行了验证。

1 双平方根偏移

为了论述需要, 首先, 简单回顾双平方根方程偏移方法。在横向变速情况下, 双平方根波场外推方程为^[8]:

$ \begin{gathered} \frac{\partial P\left(t, \mathit{\boldsymbol{X}}_{\mathrm{s}}, \mathit{\boldsymbol{X}}_{\mathrm{g}}, z\right)}{\partial z}=\left\{\left[\frac{1}{v^2\left(\mathit{\boldsymbol{X}}_{\mathrm{s}}, z\right)}-\left(\frac{\partial t}{\partial \mathit{\boldsymbol{X}}_{\mathrm{s}}}\right)^2\right]^{\frac{1}{2}}+\right. \\ \left.\left[\frac{1}{v^2\left(\mathit{\boldsymbol{X}}_{\mathrm{g}}, z\right)}-\left(\frac{\partial t}{\partial \mathit{\boldsymbol{X}}_{\mathrm{g}}}\right)^2\right]^{\frac{1}{2}}\right\} \frac{\partial P\left(t, \mathit{\boldsymbol{X}}_{\mathrm{s}}, \mathit{\boldsymbol{X}}_{\mathrm{g}}, z\right)}{\partial t} \end{gathered} $

(1)

式中: X_s=(x_s, y_s), X_g=(x_g, y_g)分别为炮点和检波点坐标; P(t, X_s, X_g, z)为z深度t时刻波场; v(X_s, z)和v(X_g, z)分别表示炮、检点层速度。ZHANG等^[3]指出该方程无成像假频的条件是需要满足成像间隔为检波点间距的一半。对于双平方根方程偏移, 直接在中点-半偏移距域实现更为自然。转化方程(1)到中点-半偏移距域, 并将波场表示为频率-波数域:

$ \begin{aligned} & \frac{\partial \widetilde{P}\left(\omega, \mathit{\boldsymbol{K}}_{\mathrm{m}}, \mathit{\boldsymbol{K}}_{\mathrm{h}} ; z\right)}{\partial z}=-\mathrm{i} k_z \widetilde{P}\left(\omega, \mathit{\boldsymbol{K}}_{\mathrm{m}}, \mathit{\boldsymbol{K}}_{\mathrm{h}} ; z\right)- \\ & \frac{\mathrm{i} k_0}{k_{z_{\mathrm{s}}}} \mathrm{FT}_{\mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}}}[\omega \Delta s(\mathit{\boldsymbol{M}}-\mathit{\boldsymbol{H}}, z) P(\omega, \mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}}, z)]- \\ & \frac{\mathrm{i} k_0}{k_{z_{\mathrm{g}}}} \mathrm{FT}_{\mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}}}[\omega \Delta s(\mathit{\boldsymbol{M}}+\mathit{\boldsymbol{H}}, z) P(\omega, \mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}}, z)] \end{aligned} $

(2)

式中: ω表示圆频率; M=(X_s+X_g)/2和H=(X_g－X_s)/2分别表示中点和半偏移距; FT表示傅里叶变换; Δs为炮点方向坐标增量; k₀=ω/v(z), k_z, K_m=(k_mx, k_my), K_h=(k_hx, k_hy)分别表示总波数、深度方向波数、中心点方向波数和半偏移距方向波数, 且有:

$ \begin{array}{l} k_{z_{\mathrm{s}}}= & \sqrt{\left(\frac{\omega}{v(z)}\right)^2-\left(\frac{\mathit{\boldsymbol{K}}_{\mathrm{m}}-\mathit{\boldsymbol{K}}_{\mathrm{h}}}{2}\right)^2} \\ k_{z_{\mathrm{g}}}= & \sqrt{\left(\frac{\omega}{v(z)}\right)^2-\left(\frac{\mathit{\boldsymbol{K}}_{\mathrm{m}}+\mathit{\boldsymbol{K}}_{\mathrm{h}}}{2}\right)^2} \\ k_z \equiv- & \operatorname{sign}(\omega)\left[\sqrt{\left(\frac{\omega}{v(z)}\right)^2-\left(\frac{\mathit{\boldsymbol{K}}_{\mathrm{m}}+\mathit{\boldsymbol{K}}_{\mathrm{h}}}{2}\right)^2}+\right. \\ & \left.\sqrt{\left(\frac{\omega}{v(z)}\right)^2-\left(\frac{\mathit{\boldsymbol{K}}_{\mathrm{m}}-\mathit{\boldsymbol{K}}_{\mathrm{h}}}{2}\right)^2}\right] \end{array} $

(3)

显然炮点和检波点方向波数满足:

$ \begin{aligned} & \left|\frac{\mathit{\boldsymbol{K}}_{\mathrm{m}}-\mathit{\boldsymbol{K}}_{\mathrm{h}}}{2}\right|=\left|\mathit{\boldsymbol{K}}_{\mathrm{s}}\right| \leqslant \frac{\pi}{\Delta s} \\ & \left|\frac{\mathit{\boldsymbol{K}}_{\mathrm{m}}+\mathit{\boldsymbol{K}}_{\mathrm{h}}}{2}\right|=\left|\mathit{\boldsymbol{K}}_{\mathrm{g}}\right| \leqslant \frac{\pi}{\Delta g} \end{aligned} $

(4)

式中: Δg为检波点方向坐标增量。

对(2)式进行分解, 得:

$ \frac{\partial \widetilde{P}\left(\omega, \mathit{\boldsymbol{K}}_{\mathrm{m}}, \mathit{\boldsymbol{K}}_{\mathrm{h}} ; z\right)}{\partial z}=-\mathrm{i} k_z \widetilde{P}\left(\omega, \mathit{\boldsymbol{K}}_{\mathrm{m}}, \mathit{\boldsymbol{K}}_{\mathrm{h}} ; z\right) $

(5)

$ \begin{gathered} \frac{\partial \widetilde{P}\left(\omega, \mathit{\boldsymbol{K}}_{\mathrm{m}}, \mathit{\boldsymbol{K}}_{\mathrm{h}} ; z\right)}{\partial z}=-\frac{\mathrm{i} k_0}{k_{z_{\mathrm{s}}}} \mathrm{FT}_{\mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}}}[\Delta s(\mathit{\boldsymbol{M}}-\mathit{\boldsymbol{H}}, z) \cdot \\ P(t, \mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}} ; z)]-\frac{i_0}{k_{z_{\mathrm{g}}}} \mathrm{FT}_{\mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}}}[\Delta s(\mathit{\boldsymbol{M}}+\mathit{\boldsymbol{H}}, z) \cdot \\ P(t, \mathit{\boldsymbol{M}}, \mathit{\boldsymbol{H}} ; z)] \end{gathered} $

(6)

(5) 式是背景介质中背景波场的相移外推方程。(6)式描述了慢度摄动引起的散射场。背景场和散射场的和构成总场。对(6)式进行不同的近似, 得到不同精度的外推方程。

(5) 式和(6)式的双平方根算子将激发点和接收点同时向下外推, 并根据零时间零偏移距成像条件:

$ P\left(t=0, \mathit{\boldsymbol{K}}_{\mathrm{m}}, H=0, z\right)=\int \mathrm{d} \omega \int \mathrm{d} \mathit{\boldsymbol{K}}_{\mathrm{h}} P\left(\omega, \mathit{\boldsymbol{K}}_{\mathrm{m}}, \mathit{\boldsymbol{K}}_{\mathrm{h}}, z\right) $

(7)

实现成像值提取, 完成叠前深度偏移^[8]。

2 数据空间假频的分析及消除

在中点-半偏移距域进行波场延拓, 中点方向的采样间隔恰好为道距的一半, 根据ZHANG等^[3]推导的双平方根方程偏移成像条件假频理论, 该间隔可以无假频地进行成像, 即它保证了(7)式成像条件无假频。(2)式、(5)式和(6)式频率域传播本身不产生假频, 即该算法无算子假频。但并不是说在该域进行波场延拓没有假频的影响, 要实现无假频的波场延拓, (2)式、(5)式和(6)式还必须满足数据$\widetilde{P}$(ω, K_m, K_h; z)是无假频的采样。在时空域, 需要满足时间采样间隔Δt、中点方向采样间隔Δm和半偏移距方向采样间隔Δh均足够小, 时间和空间方向上没有假频。时间采样间隔Δt满足时间采样定理即可, 在此不做讨论。在半偏移距方向, 根据采样定理, 无空间假频需要满足:

$ f_{\max } \leqslant \frac{1}{2 \Delta t}=\frac{1}{2 \frac{\partial t}{\partial h} \Delta h} $

(8)

式中: f_max为数据最高频率。即:

$ \Delta h \leqslant \frac{1}{2 p_{\mathrm{h}} f_{\max }} $

(9)

式中: p_h=∂t/∂h为半偏移距方向的射线参数。同理, 中点方向采样间隔Δm需要满足:

$ \Delta m \leqslant \frac{1}{2 p_{\mathrm{m}} f_{\max }} $

(10)

式中: p_m为中点方向的射线参数。

(9) 式和(10)式限定了满足中点-半偏移距域双平方根延拓的空间采样间隔。对于横向缓变速介质, 在CMP道集中, 中点-半偏移距域时距关系近似满足:

$ t \approx \sqrt{t_0^2+\left(\frac{2 h}{v_{\mathrm{rms}}}\right)^2} $

(11)

式中: v_rms为均方根速度; t₀为自激自收时间; t为地震波旅行时间。将(11)式对半偏移距h求偏导数并代入(9)式, 得:

$ \Delta h \leqslant \frac{t v_{\mathrm{rms}}^2}{8 h f_{\max }} $

(12)

(12) 式为半偏移距方向无假频必须满足的条件, 即时间越小、均方根速度越小、频率越高、偏移距越大, 要求采样间距越小。即, 浅层、低速区域和大偏移距数据容易产生假频, 偏移距间距越大越容易产生假频。

在二维共偏移距道集或者三维OVT道集中, 反射数据中点方向射线参数满足:

$ p_{\mathrm{m}}=\frac{\partial t}{\partial m}=\sin \alpha $

(13)

式中: α为反射界面的时间域倾角。通常, 倾角α较小, 因此, (13)式较容易满足。在中点-半偏移距域波场延拓中, 中点间距为道距一半, 因此, (13)式的采样条件几乎总是成立的。分析当前主流的观测数据CMP道集, 在中点间距为道距一半的情况下, 对于二维规则观测数据, 半偏移距间距恰好等于炮间距, 它通常是道间距的两倍或者三倍; 对于三维正交规则观测数据, Inline方向和Crossline方向的半偏移距间距分别为炮线间距和检波线间距, 通常为数百米的量级。即现行主流观测系统, 二维数据的炮间距和三维数据的炮线间距/检波线间距都较大, 对于中高频数据非常容易违反公式(12)的采样条件, 从而产生假频。

消除数据空间假频通常有两种方法, 方法一为压制数据的高频信息, 从而使数据满足无假频采样条件; 方法二为进行数据反假频插值, 加密数据的空间采样, 从而压制假频。方法一需要采取局部滤波措施, 避免全局滤波器滤波造成波场的过度低频化, 本文采用方法一进行高频压制。方法二会加大数据量, 从而增加算法的计算量和内存需求量, 影响算法的效率。

LUMLEY等^[5]根据三角平滑滤波器的局部低通性质通过选择恰当的滤波器长度进行低通平滑滤波, 压制高频数据, 实现无假频的Kirchhoff偏移。结合三角平滑滤波器性质和(12)式CMP域半偏移距采样条件, 设计合适的滤波器, 实现双平方根方程偏移的反假频处理。根据三角平滑滤波器性质, 滤波器长度N需要满足:

$ N \geqslant \max \left(4 \frac{\partial t}{\partial h} \cdot \frac{\Delta h}{\Delta t}-1, 1\right)=\max \left(\frac{16 h \Delta h}{t v_{\mathrm{rms}}^2 \Delta t}-1, 1\right) $

(14)

该滤波器在对地震道数据进行叠加的基础上, 通过三点运算完成, 计算效率非常高^[5]。在双平方根方程偏移波场延拓之前根据公式(14)设计的滤波器对CMP输入数据进行三角平滑滤波, 压制中、浅层低速区域数据假频, 从而实现无假频波场外推。

图 1为Marmousi模型某CMP道集反假频前后对比结果, 图 1a的半偏移距间距为50 m, 可以看到浅层大偏移距位置处明显存在假频。图 1b为根据(14)式设计的三角平滑滤波算子长度滤波后的数据。浅层大偏移距数据高频得到压制, 小偏移距和中深层没有假频的同相轴都得到了很好的保持。

3 模型数据及实际资料测试 3.1 Marmousi模型数据

Marmousi模型数据共计240炮,记录长度为3 s,时间采样间隔为4 ms。单边接收,每炮96道接收,炮间距为50 m,道间距为25 m,最大偏移距为2 575 m。成像中心点间距为12.5 m,半偏移距间距为50 m。对于该模型数据,最高频率为40 Hz,浅层0.5 s均方根速度为1 750 m/s。对于该时间,根据(12)式,半偏移距需要满足Δh≤4 785/h,即对于半偏移距为250.0 m的地震道,需要满足半偏移距间距小于19.1 m。以半偏移距50 m为间距进行成像对浅层显然会出现明显假频。对于中深层2.1 s均方根速度约为2 700 m/s,需要满足Δh≤47 840/h,即对于最大半偏移距1 287.5 m的地震道,半偏移距需要小于37.0 m,对于最大半偏移距1 000 m的地震道,半偏移距间距小于47.8 m即可。

图 2a为平滑后的Marmousi速度模型, 图 2b为利用成像点间距50 m且没有进行反假频处理双平方根方程偏移剖面, 图 2c为利用本文方法进行反假频处理后的双平方根方程偏移剖面。可以看出, 图 2b存在明显的假频噪声干扰, 尤其是中浅层, 假频非常明显, 随着深度增加, 假频逐步减少。图 2c中已经基本看不到假频噪声的影响。图 2d为将偏移距间距加密到12.5 m后的成像结果, 可以看到, 加密采样后也能够消除假频的影响, 改善成像质量。图 2e至图 2g分别为图 2b至图 2d右上角放大后的显示, 可以更明显地看到假频噪声对成像效果的影响和反假频后成像效果的改善。

根据(8)式和图 2成像结果的对比可知, 本文方法和加密采样均能够很好地实现假频噪声的压制, 但是加密空间采样会加大计算内存需求和计算量, 减弱了双平方根成像方法计算高效的优势。表 1为图 2c和图 2d计算内存需求和计算时间统计。内存需求方面, 常规加密采样方法将偏移距加密为原来的4倍, 整个数据规模也增加为原来的4倍, 需求的内存随之膨胀为原来的4倍。计算量方面, 随着内存需求的增加, 计算量也增大为原来的4倍。

结合表 1计算需求数据以及图 2成像结果对比可知, 加密采样方法和本文方法在成像上均可以有效压制假频, 但在存储和计算时间上加密采样方法均随着加密采样而数据规模显著增加, 本文方法是兼顾效率和效果的方法。

3.2 实际资料

某三维地震数据Inline方向以及Crossline方向半偏移距间隔250 m, 道间距为50 m, 中心点间距为25 m, 数据最高频率60 Hz。分别利用本文反假频方法处理前、后的双平方根方程偏移成像的结果如图 3和图 4所示, 其中, 图 3a和图 4a为未进行反假频处理的双平方根方程偏移Crossline方向和Inline方向剖面, 图 3b和图 4b为采用本文方法进行反假频处理后的Crossline方向和Inline方向成像剖面。对比可以看出, 反假频处理之前, 成像剖面中浅层假频噪声严重, 有效反射波同相轴被假频覆盖。利用本文方法进行反假频压制, 在剖面1 km上下范围内的假频噪声得到了有效压制, 有效信号得到凸显。对于该数据, 常规加密采样方法需要加密半偏移距间距到25 m, 数据规模是原来的100倍, 计算内存和计算量也膨胀100倍, 已经远远超过了测试计算机的内存, 因此加密采样对于三维数据不具备可行性。

4 讨论

本文的压制假频方法以地下缓变介质的时距曲线为基础, 对于绝大多数介质是适用的。对于复杂介质, 理论上只有准确把握同相轴的曲率才能设计准确的反假频公式, 在实际处理中既不可能做到, 也没有实际意义。在实际应用过程中, 可以在假频公式中加入调节系数, 根据地下构造倾角反假频滤波器长度进行调节, 实现不同倾角构造的假频压制, 提高方法对复杂构造的适应性。

本文方法的优势在于能够很好地平衡假频的压制和计算需求, 保障偏移的效率。众所周知, 对于假频的压制, 最直接的方式是通过反假频插值加密空间采样密度。但是加密空间采样带来的问题是计算需求的急剧提高。对于三维双平方根方程偏移, 需要同时在中点-半偏移距域进行波场延拓, 这里的中点方向和半偏移距方向分别为一个两维空间, 两个方向组成一个四维数组进行波场延拓。这个四维数组需要较大的计算机内存, 如果通过反假频插值加密半偏移距方向的采样, 数据量会急剧膨胀, 造成偏移算法无法在当前主流的计算机上进行运算或者计算效率低下。结合观测系统实际情况, 适度地加密采样并利用本文反假频方法, 减小对浅层大偏移距高频数据的压制, 得到计算量和计算效果均衡的结果。

5 结论

偏移假频通常分为数据假频、算子假频和成像条件假频, 对于中点-半偏移距域双平方根方程偏移本身不会产生传播算子及成像条件假频, 数据空间假频是该成像方法假频的唯一来源。对于地震数据, 偏移距方向的射线参数与时距曲线相关, 中点方向的射线参数与界面时间域倾角相关, 当前主流的地震观测系统, 半偏移距采样间距通常是中点间距的数倍, 中点方向不容易产生假频, 主要假频来源为半偏移距方向。

本文提出了针对中点-半偏移距域双平方根方程偏移数据空间假频消除方法。该方法在偏移成像之前对地震数据进行三角平滑滤波器滤波, 压制偏移距方向的数据假频, 从而避免波场延拓过程中偏移距方向的假频输入, 实现无假频传播及成像, 提高成像品质。该方法计算效率非常高, 与偏移成像相比, 计算量几乎可以忽略不计。模型数据和实际资料测试表明, 利用本文方法进行反假频处理后的成像剖面假频噪声得到了有效压制, 明显改善了地震数据的成像质量。